% load directory:
% close all;
clearvars
load_dir = uigetdir;

warning('off','all')
load([load_dir '\Parameters.mat'])
load([load_dir filesep 'MatchedValues.mat'])
cd(load_dir)
temp=fields(Parameters);
for i = 1:length(temp)
   eval([temp{i} '= Parameters.' temp{i} ';'])
    
end
warning('on','all')
%Exp densities
x = [49.56 139.65 229.73 319.82 409.91 500.00 590.09 680.18 770.27 860.35 950.44];
x = x./1000;
rho{1} = [0.001330 0.001127 0.001170 0.001086 0.001025 0.001044 0.001092 0.001109 0.001073 0.001131 0.001318];
rho{2} = [0.001496 0.001167 0.001036 0.001028 0.001023 0.00109 0.001085 0.001086 0.00108 0.001132 0.001355];
rho{3} = [0.001448 0.001166 0.001107 0.001046 0.001079 0.001076 0.00105 0.001036 0.001094 0.001117 0.001352];
rho{4} = [0.001404 0.001232 0.001182 0.001191 0.001149 0.001123 0.001025 0.000947 0.000959 0.001095 0.001275];
rho{5} = [0.000896 0.000834 0.000972 0.001045 0.001142 0.001102 0.001025 0.001064 0.001213 0.001477 0.001859];

%Exp Viscosity dist
Pos_Vis = [45.3 129.9 214.47 299.08 383.68 468.3 552.9 637.5 722.1 806.7 891.3]./1000;
Viscosity{1} = ones(1,length(Pos_Vis));
Viscosity{2} = [1.42 1.37 1.38 1.33 1.27 1.18 1.10 1.08 1.07 1.05 1.02];
Viscosity{3} = [2.26 2.17 2.09 1.96 1.81 1.60 1.41 1.24 1.16 1.11 1.10];
Viscosity{4} = [3.31 3.19 2.99 2.69 2.51 2.20 1.86 1.58 1.39 1.23 1.19];
Viscosity{5} = [5.57 5.47 5.29 4.75 4.17 3.58 3.01 2.28 1.75 1.45 1.26];

%Exp velocity profiles
% velocity dist %
V_exp{1} = [89.36 95.25 91.04 92.53 94.96 91.71 89.52 91.19 93.27 94.89 90.34];
V_exp{2} = [59.57 70.75 72.94 72.90 73.80 73.43 76.89 78.30 79.76 81.71 75.08];
V_exp{3} = [42.88 49.75 52.10 55.21 58.20 62.43 69.45 75.81 78.90 82.67 76.96];
V_exp{4} = [26.93 29.52 32.19 35.74 40.34 47.51 59.06 75.05 85.06 91.58 90.50];
V_exp{5} = [19.48 20.88 21.55 24.19 28.51 35.54 47.49 62.55 77.77 87.25 89.18];
x_vel = [49.56 139.65 229.73 319.82 409.91 500.00 590.09 680.18 770.27 860.35 950.44];
x_vel = x_vel./1000;

V_exp_std{1} = [13.08 10.90 13.62 11.73 10.53 12.14 12.86 13.56 11.20 9.98 10.23];
V_exp_std{2} = [0.84 1.50 2.24 1.74 1.99 1.55 0.93 1.48 1.41 0.36 1.79];
V_exp_std{3} = [3.66 4.85 4.84 5.31 5.85 4.46 3.26 2.27 3.60 1.59 3.73];
V_exp_std{4} = [1.47 2.12 1.98 2.01 2.05 1.63 3.13 5.65 7.83 9.9 9.88];
V_exp_std{5} = [0.90 1.06 1.04 1.43 2.22 2.08 2.92 3.76 4.00 3.60 3.26];

V_exp_temp = cellfun(@(z) flip(z), V_exp,'un',0);
V_exp_sigm={};
for i = 1:length(V_exp_temp)
    if i == 1
        V_exp_sigm{i} = @(z) 1+0.*z;
        slope_sigm(1) = 0;
        continue
    end
    param = sigm_fit(x_vel(2:end-1),V_exp_temp{i}(2:end-1));
    temp = @(x) (param(1)+(param(2)-param(1))./(1+10.^((param(3)-x)*param(4))));
    temp = @(x) temp(x)/temp(0);
    V_exp_sigm{i} = temp;
        temp = max(abs(diff(V_exp_sigm{i}(0:0.01:1))./0.01));
    slope_sigm(i) = temp;
end


color_mark={'b','--b','c','--c','r','--r','g','--g','m','--m',};
figure; hold on;
for i = 1:length(rho); plot(x,rho{i}./trapz(x,rho{i}),color_mark{2*i-1}); plot(Pos_Vis, Viscosity{i}./trapz(Pos_Vis,Viscosity{i}), color_mark{2*i});end
legend({'control den', 'control vis','1.5 cp den', '1.5 cp vis', '2 cp den',...
    '2 cp vis', '5 cp den', '5 cp vis', '9 cp den','9 cp vis'},'location','eastoutside');
title('Cell Density and Viscosity (from exp)')

%% Matched Density (single plot)
figure; hold on; box on;
x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x = x./1000;
colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]}; count = 1;
for i = 1:5
    load([load_dir,'\tau',num2str(Tau_match(i)),'Vr',num2str(Vr_match(i)),...
    'Vl',num2str(V_0),'Drot',num2str(D_rot_dim(i)),'.mat']);
    
    [counts, centers] = hist(XX(:),11);
    counts = flip(counts);
    centers = centers * 1000;
    counts = counts./trapz(centers,counts);
    sim_den{i} = counts; sim_x{i} = centers;
    plot(centers,counts,'color',colors{count},'linewidth',1)
    plot(x*1000, rho{i}./trapz(x*1000,rho{i}),'color', colors{count}, 'linestyle', '--','linewidth',1)
    count = count + 1;
end
ylim([0.8 1.8]*10^(-3))
yticks([0.8 1 1.2 1.4 1.6 1.8]*10^(-3))
xticks([0 250 500 750 1000])
% xticklabels({''});
% yticklabels({''});
xlabel('position')
ylabel('PDF')
% legend({'control sim','control exp', '1.5 cp sim', '1.5 cp exp','2 cp sim',...
%     '2 cp exp', '5 cp sim', '5 cp exp', '9 cp sim', '9 cp exp'},...
%    'location', 'north')
cd(load_dir)
set(gca,'linewidth', 1)
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [4 4]);
cd(load_dir)
% print('Matched_Densities','-dpdf','-r1200', '-bestfit')
%% Matched Densities (subplots)
figure;
x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x = x./1000;
colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]}; count = 1;
ha = tight_subplot(1,5,0.01,0.01,0.06);
for i = 1:5
    axes(ha(i));
    %   subplot(5,1,i)
    hold on; box on;
    load([load_dir,'\tau',num2str(Tau_match(i)),'Vr',num2str(Vr_match(i)),...
        'Vl',num2str(V_0),'Drot',num2str(D_rot_dim(i)),'.mat']);
    
    [counts, centers] = hist(XX(:),11);
    counts = flip(counts);
    centers = centers * 1000;
    counts = counts./trapz(centers,counts);
    plot(centers,counts,'color',colors{count},'linewidth',1)
    plot(x*1000, rho{i}./trapz(x*1000,rho{i}),'color', colors{count}, 'linestyle', '--','linewidth',1)
    count = count + 1;
    
    ylim([0.8 1.8]*10^(-3))
    xlim([0 1000])

    if i == 1
        ylabel('PDF')
        yticks([0.8 1 1.2 1.4 1.6 1.8]*10^(-3))
        yticklabels({'0.8','1','1.2','1.4','1.6','1.8'})
    else
        yticks([])
    end
    xticks([0 250 500 750])   
    xticklabels({'0','250','500','750'})
    xlabel('position')
    set(gca,'linewidth', 1)
    daspect([10^6 1 1])
    
end


% 
% set(gcf, 'PaperUnits', 'inches');
% set(gcf, 'PaperSize', [8 4]);
cd(load_dir)
% print('Matched_Densities','-dpdf','-r1200', '-bestfit')
%%
figure; hold on; box on;
x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x_sim = x./1000;
colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]}; count = 1;
for i = 1:5
    load([load_dir,'\tau',num2str(Tau_match(i)),'Vr',num2str(Vr_match(i)),...
    'Vl100','Drot',num2str(D_rot_dim(i)),'.mat']);
    
    [counts, centers] = hist(XX(:),11);
    counts = flip(counts);
    centers = centers * 1000;
    counts = counts./trapz(centers,counts);
    
    sim_v = linspace(Vr_match(i), 100, length(x));
    sim_v = sim_v./trapz(x, sim_v);
    exp_v = V_exp{i}./trapz(x, V_exp{i});
    
    plot(centers,(counts.*sim_v)./trapz(centers,(counts.*sim_v)),'color',colors{count},'linewidth',2)
    plot(x, ((rho{i}./trapz(x,rho{i})).*exp_v)./trapz(x, ((rho{i}./trapz(x,rho{i})).*exp_v)),'color', colors{count}, 'linestyle', '--','linewidth',2)
    count = count + 1;

end
xticks([0 250 500 750 1000])
xlabel('position')
ylabel('schnitzer departure')
% legend({'control sim','control exp', '1.5 cp sim', '1.5 cp exp','2 cp sim',...
%    '2 cp exp', '5 cp sim', '5 cp exp', '9 cp sim', '9 cp exp'},...
%    'location', 'north')
cd(load_dir)
% cd('C:\Users\dwalka01\Desktop\Viscotaxis_Simulations\Figures\Figure3\matlab')
% print('SchnitzerDeparture','-dpdf','-r1200')

%%
load('C:\Users\dwalka01\Desktop\Viscotaxis_Simulations\CorrectedExpTracks\Turning_Rates(33-66).mat')
row = find(V_r(:,1) == Vr_match(5));
cmap = cool(256);
figure; hold on; box on;
colormap cool

ind = [1:length(V_r(:,1))];

for ii = ind
   load([load_dir,'\tau',num2str(Tau_visc(row,ii)),'Vr',num2str(V_r(row,ii)),...
    'Vl',num2str(V_0),'Drot',num2str(D_rot_dim(row,ii)),'.mat'], 'XX'); 
    [counts, centers] = hist(XX(:),11);counts = counts./trapz(centers,counts);
    temp = log10(1./Tau_visc(1,:));temp = temp + abs(min(temp)); temp = temp./max(temp);

    plot(centers, flip(counts), 'color', cmap(floor(255*temp(ii))+1,:))
        
    if ii == ind(end)
        V = (1+centers*(V_r(row,ii)-100)/100);
        k = 1./trapz(centers, 1./V);
        rho_s = k./V;
        plot(centers, flip(rho_s), '--k', 'linew', 2)
    end
    
end
cbar = colorbar;
caxis([log10(1./Tau_visc(1,end)) log10(1./Tau_visc(1,1))])
cbar.Ticks = [-3 -2 -1];
xticks([0 0.25 0.5 0.75 1])
xlabel('position')
ylabel('Density (PDF)')
title(cbar, 'log10 Turning rate (1/\tau)')
% xticklabels({''});
% yticklabels({''});
% xlabel('position')
% ylabel('PDF')
% legend({'control sim','control exp', '1.5 cp sim', '1.5 cp exp','2 cp sim',...
%     '2 cp exp', '5 cp sim', '5 cp exp', '9 cp sim', '9 cp exp'},...
%    'location', 'north')
cd(load_dir)
set(gca,'linewidth', 1.5)
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [6 4]);

cd(load_dir)
%print('phasematrixverticallinecut', '-dpdf', '-r1200', '-bestfit')
%% Domains (viscophobic viscophilic)
close all
PM = zeros(size(B));
PM_x = PM; PM_V = PM; PM_V_norm = PM;
x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x_sim = x./1000;
V_r = Parameters.V_r;
for row=1:length(B(:,1))
    
for ii = 1:length(B(row,:))
% load([load_dir,'\B',num2str(B(ii)),'Vr',num2str(V_r(ii)),'Vl',num2str(V(ii)),'.mat']);
load([load_dir,'\tau',num2str(Tau_visc(row,ii)),'Vr',num2str(V_r(row,ii)),...
    'Vl',num2str(V_0),'Drot',num2str(D_rot_dim(row,ii)),'.mat']);
temp=XX(:);

[counts, centers] = hist(XX(:),11);counts = counts./trapz(centers,counts);
V = (1+centers*(V_r(row,ii)-100)/100);
k = 1./trapz(centers, 1./V);
rho_s = k./V;
% PM(row, ii) = trapz(centers,counts.*log(counts./rho_s));
PM(row, ii) = sqrt(mean((counts./trapz(centers,counts)-rho_s).^2));
% 

% temp=temp(:);
% ind=find(temp>0.5);
% PM(row,ii)=(2*length(ind)-length(temp))/(length(temp));
% PM(ii) = mean(XX(:));
end
row
end
PM=PM';
PM=flip(PM,2);

% PM_x=PM_x'; PM_x=flip(PM_x,2);
% PM_V=PM_V'; PM_V=flip(PM_V,2);
% PM_V_norm=PM_V_norm'; PM_V_norm=flip(PM_V_norm,2);

% save([load_dir filesep 'PhaseMatrix.mat'],'PM')
% save([load_dir filesep 'PhaseMatrix_schnitzer.mat'],'PM')
% save([load_dir filesep 'SkewnessMatrix.mat'],'PM')
% save([load_dir filesep 'Int_xrhov_Matrix.mat'],'PM')
% save([load_dir filesep 'x_V_Vnorm_Skewnessmatrices.mat'],'PM_x', 'PM_V', 'PM_V_norm')

% save([load_dir filesep 'KullbackDivergence(rho.rho_s).mat'],'PM')
save([load_dir filesep 'RMS_diff(rho.rho_s).mat'],'PM')

%% Log Plotting

load([load_dir filesep 'RMS_diff(rho.rho_s).mat'],'PM');

load('C:\Users\dwalka01\Desktop\Viscotaxis_Simulations\CorrectedExpTracks\Turning_Rates(33-66).mat')
% grad_eta = [0, 0.79, 2.2, 3.4, 7.2];
for i = 1:length(Viscosity)
   temp = interp1(Pos_Vis, Viscosity{i}, [0.33 0.66]);
   grad_eta_third(i)= (temp(1)-temp(2))./0.333;
   grad_eta_full(i) = Viscosity{i}(1) - Viscosity{i}(end);
   grad_eta(i) = temp(1)-temp(2);
end
scaling_arg = nanmean(grad_eta_third./grad_eta_full);

clim = 0.3;

x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x = x./1000;
accum_exp = cellfun(@(z) (2*trapz(x(1:6),z(1:6))-trapz(x,z))./trapz(x,z), rho, 'un', 0);
accum_exp = [accum_exp{:}];

% h_filter = fspecial('gaussian',5,2);
% PM = imfilter(PM, h_filter, 'replicate')

figure('color',[1 1 1]); 
imagesc([(list_V_r(1)-list_V_r(end))./list_V_r(1) (list_V_r(1)-list_V_r(1))./list_V_r(1)],log10([1/list_Tau_visc(1) 1/list_Tau_visc(end)]),PM); hold on

% M1 = contour((list_V_r-list_V_r(end))./list_V_r(1),log10(1./list_Tau_visc),PM, 30,'k');

ind = linspace(0,1,256);
color_formap = ([100 ; 100 ; 200]./256).*0.8;
cmap = color_formap + ind.*([1; 1 ;1]-color_formap);
cmap = cmap';
cmap = flip(cmap);
colormap(cmap);
cbar = colorbar;
cbar.LineWidth = 1;
cbar.Color = 'k';
% 

% cmap = colormap('redblue');
% colormap(cmap);
% colorbar
% hold on; 

colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]};
% Tau_visc_exp = [.0015 .0154 0.0413 0.0429 0.0668];
% Tau_visc_exp = [0.0011371343134154, 0.0238302174134351, 0.0483465830594653, 0.0501471376164769, 0.0688445939937381];
Tau_visc_exp = mean_omega;

D_rot_exp = 0.15;
% Tau_visc_error = [0.0058 0.0059 0.0062 0.0082 0.0053];
% Tau_visc_error = [0.00307990740137862, 0.00354201504964154, 0.00381563714908374, 0.00518531019183087, 0.00345658234189232];
Tau_visc_error = std_omega;

V_r_error_temp = cellfun(@(x) sqrt(x(end).^2 + x(1).^2), V_exp_std, 'un', 0); V_r_error_temp = [V_r_error_temp{:}];

Tau_visc_min = Tau_visc_exp - Tau_visc_error; Tau_visc_max = Tau_visc_exp + Tau_visc_error;
ind = find(Tau_visc_min < 0); Tau_visc_min(ind) = 0.00001;
% 
% plot(Contour_03(1,:), Contour_03(2,:), '--k')
% plot(Contour_06(1,:), Contour_06(2,:), 'k')
% plot(Contour_12(1,:), Contour_12(2,:), '--k')

% hline(log10(0.06), '--k')

exponent = 0.934;

eta_over_eta_0 = 1.01:0.01:20;
eta_over_eta_0_exp = cellfun(@(x) x(1)/x(end), Viscosity, 'un',0);
eta_over_eta_0_exp = [eta_over_eta_0_exp{:}];
W=1000;
alpha = 33;
y = scaling_arg*(2*alpha)/W*(1-(2./(eta_over_eta_0+1)));
y_exp = scaling_arg*(2*alpha)/W*(1-(2./(eta_over_eta_0_exp+1)));
x = 1-1./eta_over_eta_0.^(exponent);
x_exp = 1-1./eta_over_eta_0_exp.^(exponent);
plot(x,log10(y),'-k', 'linew',2);
for i =1:length(x_exp);plot(x_exp(i),log10(y_exp(i)),'dk','markerfacecolor',colors{i});end

 yticks([-3 -2 -1 -0.06]);
 yticklabels({'10^{-3}','10^{-2}','10^{-1}', '10^{0}'})


V_rr = [];
for i = 1:length(Tau_visc_min)
    V_rr(i) = (V_exp_sigm{i}(0)-V_exp_sigm{i}(1))./V_exp_sigm{i}(0);
    V_r_temp = (V_exp{i}(end)-V_exp{i}(1))./V_exp{i}(end);
    V_r_error(i) = sqrt((V_r_error_temp(i)./(V_exp{i}(end)-V_exp{i}(1))).^2 + (V_exp_std{i}(end)./V_exp{i}(end)).^2).*V_r_temp;
    x1 = ones(1,100)*V_rr(i);
    y1 = linspace(log10(Tau_visc_min(i)),log10(Tau_visc_max(i)),100);
    y2 = ones(1,100)*log10(Tau_visc_exp(i));
    x2 = linspace(V_rr(i)-V_r_error(i),V_rr(i)+V_r_error(i),100);
    line(x1, y1, 'color', colors{i}, 'linewidth', 2)
    line(x2, y2, 'color', colors{i}, 'linewidth', 2)
end
for i=1:5
% plot(V_rr(i),log10(Tau_visc_exp(i)),'o','markersize',6,'linew',2,'color',colors{i}, 'markerfacecolor', cmap(floor(256*accum_exp(i)/(2*clim))+128,:));
plot(V_rr(i),log10(Tau_visc_exp(i)),'o','markersize',6,'linew',2,'color',colors{i}, 'markerfacecolor', 'w');

Vr_match_new =  (100 - Vr_match_exact(i))./100;% (100 - [100, 81.3793, 56.5517, 31.7241, 22.41])./100;
Tau_visc_match = 1./Tau_match(i); % 1./[1000 398.11 31.62 31.62 10];
% plot(Vr_match_new,log10(Tau_visc_match),'sk','markersize',8,'linew',2,'color',colors{i},'markerfacecolor',cmap(floor(256*accum_exp(i)/(2*clim))+128,:));
plot(Vr_match_new,log10(Tau_visc_match),'sk','markersize',8,'linew',2,'color',colors{i},'markerfacecolor','w');

end 


% yticks([log10(0.006) log10(0.06)]);
% yticklabels({'10^{-1}','10^0'})
xlabel('$\Delta v/v_\mathrm{max}$','Interpreter','Latex','Fontsize',16);
ylabel('$1/\tau_\mathrm{visc}$','Interpreter','Latex','Fontsize',16);
% ylabel('$\textrm{Pe}_r$','Interpreter','Latex','Fontsize',16);

% caxis([-clim clim])
% caxis([-3 -0.657])
% caxis([0 0.8])

% title('viscophobic (Blue) and viscophilic (Red) domains')
% legend(plt, {'experimental','matched'},'fontsize',16)
set(gca,'Ydir','normal')
% daspect([1 5 1])
% cd('C:\Users\dwalka01\Desktop\Viscotaxis_Simulations\Figures\Figure3\matlab')
cd(load_dir)
daspect([1 2.9 1])
% print('PhaseMap_accum(sidebar).pdf','-dpdf','-r1200')
% print('PhaseMap_accum.pdf','-dpdf','-r1200')
%  print('PhaseMap_schnitzer.pdf','-dpdf','-r1200')
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1;
 title("rho|rho_s")
     print('PhaseMap_RMS(rho.rho_s).pdf','-dpdf','-r1200')


%% Lin Plotting

% load([load_dir filesep 'PhaseMatrix_schnitzer.mat'],'PM')
% PM_schnitz = PM;
load([load_dir filesep 'PhaseMatrix.mat'],'PM')
% load([load_dir filesep 'Contours.mat'])
clim = 0.3;


x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x = x./1000;
accum_exp = cellfun(@(z) (2*trapz(x(1:6),z(1:6))-trapz(x,z))./trapz(x,z), rho, 'un', 0);
accum_exp = [accum_exp{:}];

% h_filter = fspecial('gaussian',5,2);
% PM = imfilter(PM, h_filter, 'replicate')

figure('color',[1 1 1]); 
imagesc([(list_V_r(1)-list_V_r(end))./list_V_r(1) (list_V_r(1)-list_V_r(1))./list_V_r(1)],[1/list_Tau_visc(1) 1/list_Tau_visc(end)],PM); hold on

% M1 = contour((list_V_r-list_V_r(end))./list_V_r(1),log10(1./list_Tau_visc),PM, [0 0],'k');
 cmap = colormap('redblue');
ind = linspace(0,1,256);
% cmap = [0; 0 ;1] + ind.*([1; 1 ;1]-[0; 0; 1]);
% cmap = cmap';
colormap(cmap);
colorbar
hold on; 

colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]};
% Tau_visc_exp = [.0015 .0154 0.0413 0.0429 0.0668];
Tau_visc_exp = [0.0011371343134154, 0.0238302174134351, 0.0483465830594653, 0.0501471376164769, 0.0688445939937381];
D_rot_exp = 0.15;
% Tau_visc_error = [0.0058 0.0059 0.0062 0.0082 0.0053];
Tau_visc_error = [0.00307990740137862, 0.00354201504964154, 0.00381563714908374, 0.00518531019183087, 0.00345658234189232];
V_r_error_temp = cellfun(@(x) sqrt(x(end).^2 + x(1).^2), V_exp_std, 'un', 0); V_r_error_temp = [V_r_error_temp{:}];

Tau_visc_min = Tau_visc_exp - Tau_visc_error; Tau_visc_max = Tau_visc_exp + Tau_visc_error;
ind = find(Tau_visc_min < 0); Tau_visc_min(ind) = 0.00001;
% 
% plot(Contour_03(1,:), Contour_03(2,:), '--k')
% plot(Contour_06(1,:), Contour_06(2,:), 'k')
% plot(Contour_12(1,:), Contour_12(2,:), '--k')

% hline(log10(0.06), '--k')
yticks([-0.1 -0.05 0 0.05 0.1]);
V_r = [];
for i = 1:length(Tau_visc_min)
    V_r(i) = (V_exp_sigm{i}(0)-V_exp_sigm{i}(1))./V_exp_sigm{i}(0);
    V_r_temp = (V_exp{i}(end)-V_exp{i}(1))./V_exp{i}(end);
    V_r_error(i) = sqrt((V_r_error_temp(i)./(V_exp{i}(end)-V_exp{i}(1))).^2 + (V_exp_std{i}(end)./V_exp{i}(end)).^2).*V_r_temp;
    x1 = ones(1,100)*V_r(i);
    y1 = linspace(Tau_visc_min(i),Tau_visc_max(i),100);
    y2 = ones(1,100)*Tau_visc_exp(i);
    x2 = linspace(V_r(i)-V_r_error(i),V_r(i)+V_r_error(i),100);
    line(x1, y1, 'color', colors{i}, 'linewidth', 2)
    line(x2, y2, 'color', colors{i}, 'linewidth', 2)
end
for i=1:5
% plot(V_r(i),log10(Tau_visc_exp(i)),'o','markersize',6,'linew',2,'color',colors{i}, 'markerfacecolor', cmap(floor(256*accum_exp(i)/(2*clim))+128,:));
plot(V_r(i),Tau_visc_exp(i),'o','markersize',6,'linew',2,'color',colors{i}, 'markerfacecolor', 'w');

Vr_match_new =  (100 - Vr_match_exact(i))./100;% (100 - [100, 81.3793, 56.5517, 31.7241, 22.41])./100;
Tau_visc_match = 1./Tau_match(i); % 1./[1000 398.11 31.62 31.62 10];
% plot(Vr_match_new,log10(Tau_visc_match),'sk','markersize',8,'linew',2,'color',colors{i},'markerfacecolor',cmap(floor(256*accum_exp(i)/(2*clim))+128,:));
plot(Vr_match_new,Tau_visc_match,'sk','markersize',8,'linew',2,'color',colors{i},'markerfacecolor','w');

end 

%  yticks([-3 -2 -1 -0.06]);
%  yticklabels({'10^{-3}','10^{-2}','10^{-1}', '10^{0}'})

%   
 
% yticks([log10(0.006) log10(0.06)]);
% yticklabels({'10^{-1}','10^0'})
xlabel('$\Delta v/v_\mathrm{max}$','Interpreter','Latex','Fontsize',16);
ylabel('$1/\tau_\mathrm{visc}$','Interpreter','Latex','Fontsize',16);
% ylabel('$\textrm{Pe}_r$','Interpreter','Latex','Fontsize',16);

caxis([-clim clim])
% caxis([-0.3 0.3])

% title('viscophobic (Blue) and viscophilic (Red) domains')
% legend(plt, {'experimental','matched'},'fontsize',16)
set(gca,'Ydir','normal')
% daspect([1 5 1])
% cd('C:\Users\dwalka01\Desktop\Viscotaxis_Simulations\Figures\Figure3\matlab')
cd(load_dir)
% daspect([1 2.9 1])
print('PhaseMap_accum(sidebar).pdf','-dpdf','-r1200')
%   print('PhaseMap_accum.pdf','-dpdf','-r1200')
%  print('PhaseMap_schnitzer.pdf','-dpdf','-r1200')











